classdef GetRFSpectrum < Table1.NMRGJobs.NMRGJob
    methods
        function obj = GetRFSpectrum(nmrg)
            obj@Table1.NMRGJobs.NMRGJob(nmrg, 'GetRFSpectrum');
            obj.defineParameter();
        end
        
        function obj = defineParameter(obj)
            obj.paramParser=inputParser;
            obj.paramParser.addParameter('initialField', 6.0, @(x) x>0);
            obj.paramParser.addParameter('resonantFrequency', 100e3, @(x) x>0);
            obj.paramParser.addParameter('drivingAmplitude', 0.1, @(x) x>0);
            obj.paramParser.addParameter('bzMin', 0.0, @(x) x>=0);
            obj.paramParser.addParameter('bzMax', 7.0, @(x) x>=0);
            obj.paramParser.addParameter('npoint', 100, @(x) x>0);
            obj.paramParser.addParameter('average', 100, @(x) x>0);
            obj.paramParser.addParameter('signalOscillator', 1, @Miscellaneous.isZIindex4);
            obj.paramParser.addParameter('order', 4, @(x) ismember(x, 1:8));
            obj.paramParser.addParameter('bandwidth', 100, @(x) x>0);
            obj.paramParser.addParameter('fittingWidth', 1.0, @(x) x>0);
            obj.paramParser.addParameter('isFixBz0', false, @islogical);
            obj.paramParser.addParameter('isPlot',   true, @islogical);
            obj.paramParser.addParameter('description', 'no.', @ischar);            

        end
        %%
        function obj = preprocess(obj)
            param = obj.dataPackage.parameters;
            
            obj.boss.setRFMagnetometer('field',            param.initialField, ...
                'frequency',        param.resonantFrequency, ...
                'amplitude',        param.drivingAmplitude, ...
                'signalOscillator', param.signalOscillator, ...
                'order',            param.order, ...
                'bandwidth',        param.bandwidth);
        end
        
        %%
        function obj = dataAcquisition(obj)
            param = obj.dataPackage.parameters;
            
            measurement_data.start = datetime('now');
            obj.boss.rbSignal.unsubscribe_signal();

            
            obj.boss.rbSignal.subscribe_signal(ziDemodTime.R_avg, [2 1 1]);
            obj.boss.rbSignal.subscribe_signal(ziDemodTime.Theta_avg, [2 1 2]);
            measurement_data.signal =  obj.boss.rbSignal.getParameterDependence(...
                obj.boss.ziSlave, ...
                ziSweepNode.AuxOut3_Offset, ...
                param.bzMin, ...
                param.bzMax, ...
                'npoint', param.npoint, ...
                'averaging_sample', param.average);
            measurement_data.finish = datetime('now');
            
            obj.dataPackage.measurement_data = measurement_data;
            
        end

        %%
        function obj = analyze(obj)
            param = obj.dataPackage.parameters;
            measuredData = obj.dataPackage.measurement_data;
            
            data = measuredData.signal.getNodeData(1);
            bzList = data.x; signalList= data.y;
            
            x = bzList;
            [maxY, pkIdx] = max(signalList);
            y2 = (signalList./maxY).^2;
            [xData, yData] = prepareCurveData( x, y2 );
            
            center0 = bzList(pkIdx);
            range = param.initialField + [-param.fittingWidth,  param.fittingWidth];
            range_width = range(2) - range(1);
            excludeIdx = find( x<range(1) | x>range(2));
            
            if param.isFixBz0
                rangeLeft = param.initialField;
                center = param.initialField;
                rangeRight = param.initialField;
            else
                rangeLeft = range(1);
                center = center0;
                rangeRight = range(2);
            end
            
            ft = fittype( 'a*g*g/((x-x0)^2+g*g)', 'independent', 'x', 'dependent', 'y' );
            excludedPoints = excludedata( xData, yData, 'Indices', excludeIdx);
            opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
            opts.Display = 'Off';
            opts.Lower =      [0.5 0.0             rangeLeft];
            opts.StartPoint = [1.0 0.5*range_width center];
            opts.Upper =      [2.0 1.0*range_width rangeRight];
            opts.Exclude = excludedPoints;
            
            [analysis_result.result, analysis_result.gof] = fit( xData, yData, ft, opts );
            errList = confint(analysis_result.result);
            
            analysis_result.bz0 = analysis_result.result.x0;
            analysis_result.bz0_err = errList(:, 3);
            analysis_result.width = analysis_result.result.g;
            analysis_result.width_err = errList(:, 2);
            analysis_result.hight = sqrt(analysis_result.result.a)*maxY;
            analysis_result.hight_err = sqrt(errList(:, 1))*maxY;
            analysis_result.range = range;
            analysis_result.amplitude = maxY;
            
            obj.dataPackage.analysis_result = analysis_result;
        end
        
        %%
        function obj = postprocess(obj)
            L = log4m.getLogger;
            
            param = obj.dataPackage.parameters;
            measuredData = obj.dataPackage.measurement_data;
            analysis = obj.dataPackage.analysis_result;
            
            % Post operations
            obj.boss.bz0.setDC(analysis.bz0); pause(0.5);
            L.info('GetRFSpectrum::postprocess', sprintf('Main field Control voltage is set to [%5.3f] V', analysis.bz0));
            phase0 = obj.boss.rbSignal.autoZeroPhase();
            L.info('GetRFSpectrum::postprocess', sprintf('Demod Phase is set to [%5.3f]', phase0));
            obj.boss.by.disable();
            
            % summary
            summary.start = datestr(measuredData.start);
            summary.main_frequency = param.resonantFrequency;
            summary.phase0 = phase0;
            summary.bz0 = analysis.bz0;
            summary.bz0_err = diff(analysis.bz0_err)/2;
            summary.width = analysis.width;
            summary.width_err = diff(analysis.width_err)/2;
            summary.hight = analysis.hight;
            summary.hight_err = diff(analysis.hight_err)/2;
            summary.rsquare = analysis.gof.rsquare;
            summary.amplitude = analysis.amplitude;
            summary.finish = datestr(measuredData.finish);
            
            obj.dataPackage.summary = summary;
            obj.boss.SystemInfo.MainField = summary;
            
            % Store Data
            getSession;
            sess.addData('GetRFSpectrum', obj.dataPackage, param.description);
        end

        %%
        function obj = presentation(obj, viewer)
            param        = obj.dataPackage.parameters;
            measuredData = obj.dataPackage.measurement_data;
            analysis     = obj.dataPackage.analysis_result;
            summary      = obj.dataPackage.summary;
            
            if nargin == 1
                f = figure('Name', sprintf('RFResonance[%s]',datestr(measuredData.start)), 'WindowStyle', 'docked');
                ax1 = subplot(2, 1, 1, 'Parent', f); ax2 = subplot(2, 1, 2, 'Parent', f);
            else
                Table1.NMRGJobs.NMRGJob.updateViewer(viewer, param, summary);            
                ax1 = subplot(2, 1, 1, 'Parent', viewer.gui.uiHandle.panelPlot);
                ax2 = subplot(2, 1, 2, 'Parent', viewer.gui.uiHandle.panelPlot);
            end
            
            amplitude = analysis.amplitude;
            range = analysis.range;
            bz0 = analysis.bz0;
            Gamma = analysis.width;
            hight = analysis.hight;
            xList = linspace(range(1), range(2), 101);
            yList = sqrt(analysis.result(xList))*amplitude;
            rangeCenter = mean(range);
            
            yPhaseList = (atan((xList-bz0)/Gamma) - atan((rangeCenter-bz0)/Gamma))/pi*180;
            
            dataR = measuredData.signal.getNodeData(1); bzList = dataR.x; rList=dataR.y;
            dataTheta = measuredData.signal.getNodeData(2); phaseList = dataTheta.y;
            
            plot(ax1, bzList, rList, 'b.-', xList, yList, 'r-'); grid(ax1, 'on');
            xlabel(ax1, 'Bz0 Voltage (V)'); ylabel(ax1, 'Amplitude (V)');
            title(ax1, sprintf('Resoance Bz0=%5.4f V, \\Gamma=%5.3fV, Vmax=%5.3e V', bz0, Gamma, hight), 'Interpreter', 'tex');
            plot(ax2, bzList, phaseList/pi*180, 'b.-', xList, yPhaseList, 'r-'); grid(ax2, 'on');
            xlabel(ax2, 'Bz0 Voltage (V)'); ylabel(ax2, 'Phase (degree)');
        end

    end
    %%
    methods(Static)
        function res = compare(objList, varargin)
            p = inputParser;
            p.addParameter('isPlot', true, @islogical);
            p.addParameter('viewer', [], @(x)  isa(x, 'Table1.GUI.DataViewer') || isempty(x));
            p.parse(varargin{:});
            
            res = compare@Table1.NMRGJobs.NMRGJob(objList, 'viewer', p.Results.viewer);

            if p.Results.isPlot
                startTime = arrayfun(@(t) t.start, res.summary, 'UniformOutput', false);
                
                if isempty(p.Results.viewer)
                    f = figure('Name', sprintf('RFResonance'));
                    ax1 = subplot(2, 1, 1, 'Parent', f); ax2 = subplot(2, 1, 2, 'Parent', f);
                else
                    ax1 = subplot(2, 1, 1, 'Parent', p.Results.viewer.gui.uiHandle.panelPlot);
                    ax2 = subplot(2, 1, 2, 'Parent', p.Results.viewer.gui.uiHandle.panelPlot);
                end
                
                dataR=arrayfun(@(d) d.signal.getNodeData(1), res.measuredData);
                bzList=cell2mat(arrayfun(@(r) r.x', dataR, 'UniformOutput', false));
                rList=cell2mat(arrayfun(@(r) r.y', dataR, 'UniformOutput', false));
                
                dataTheta=arrayfun(@(d) d.signal.getNodeData(2), res.measuredData);
                phaseList=cell2mat(arrayfun(@(t) t.y', dataTheta, 'UniformOutput', false));
                
                plot(ax1, bzList, rList, '.-'); grid(ax1, 'on');
                xlabel(ax1, 'Bz0 Voltage (V)'); ylabel(ax1, 'Amplitude (V)');
                title(ax1, 'RF Resonance', 'Interpreter', 'tex');
                legend(ax1, startTime, 'Location', 'eastoutside');
                plot(ax2, bzList, phaseList/pi*180, '.-'); grid(ax2, 'on');
                xlabel(ax2, 'Bz0 Voltage (V)'); ylabel(ax2, 'Phase (degree)');
                legend(ax2, startTime, 'Location', 'eastoutside');
            end
        end
    end
end


